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We present an efficient learning algorithm for the problem of training neural networks 
with discrete synapses, a well-known hard (NP-complete) discrete optimization problem. 
The algorithm is a variant of the so-called Max-Sum (MS) algorithm. In particular, we 
show how, for bounded integer weights with q distinct states and independent concave a 
priori distribution (e.g. h regularization), the algorithm’s time complexity can be made to 
scale as O {N log N) per node update, thus putting it on par with alternative schemes, such 
as Belief Propagation (BP), without resorting to approximations. Two special cases are of 
particular interest: binary synapses W € { — 1,1} and ternary synapses W € { — 1, 0,1} with 
Iq regularization. The algorithm we present performs as well as BP on binary perceptron 
learning problems, and may be better suited to address the problem on fully-connected two- 
layer networks, since inherent symmetries in two layer networks are naturally broken using 
the MS approach. 
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I. INTRODUCTION 

The problem of training an artificial, feed-forward neural network in a supervised way is a well-known 
optimization problem, with many applications in machine learning, inference etc. In general terms, the 
problem consists in obtaining an assignment of “synaptic weights” (i.e. the parameters of the model) such 
that the device realizes a transfer function which achieves the smallest possible error rate when tested on 
a given dataset of input-output examples. Time is usually assumed to be discretized. In a single-layer 
network, the transfer function is typically some non-linear function (e.g. a sigmoid or a step function) of 
the scalar product between a vector of inputs and the vector of synaptic weights. In multi-layer networks, 
many single-layer units operate in parallel on the same inputs, and their outputs provide the input to 
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other similar (with a varying degree of similarity) units, until the last layer is reached. 

The most popular and successful approaches to these kind of optimization problems are typically 
variants of the gradient descent algorithm, and in particular the back-propagation algorithm [1]. On 
single-layer networks with simple non-linearities in their output functions these algorithms can even be 
shown to achieve optimal results in linear time [2]; on multi-layer networks these algorithms suffer from 
the usual drawbacks of gradient descent (mostly the presence of local minima, and slow convergence under 
some circumstances). 

On the other hand, gradient descent can only be applied to continuous problems. If the synaptic 
weights are restricted to take only discrete values, the abovementioned family of methods can not be 
applied; in fact, it is known that even the simplest version of the problem (classification using a single¬ 
layer network) becomes computationally hard (NP-complete) in the worst-case scenario [310]. However, 
some theoretical properties of the networks, such as the storage capacity (i.e. the amount of information 
which can be effectively stored in the device by setting the synaptic weights), are only slightly worse in 
the case of discrete synapses, and other properties (e.g. robustness to noise and simplicity) would make 
them an attractive model for practical applications. Indeed, some experimental results ISHZ!, as well as 
arguments from theoretical studies and computer simulations [EHII], suggest that long term information 
storage may be achieved by using discrete — rather than continuous — synaptic states in biological neural 
networks. 

Therefore, the study of neural network models with discrete weights is interesting both as a hard 
combinatorial optimization problem and for its potential applications, in practical implementations as 
well as for modeling biological networks. On the theoretical side, some light has been shedded upon the 
origin of the computational hardness in these kind of problems by the study of the space of the solutions 
by means of methods derived from Statistical Physics approaches nans]: in brief, most solutions are 
isolated, i.e. far from each other, and the energy landscape is riddled with local minima which tend to 
trap purely local search methods, which thus show very poor performance. On the application side, a 
family of heuristic algorithms, derived from the cavity method, have been devised, which exhibit very 
good performance on random instances, both in terms of solution time and in terms of scaling with the 
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size of the problem. 

In particular, it was first shown in |T3j that a version of the Belief Propagation (BP) algorithm [TB] 
with the addition of a reinforcement term was able to efficiently solve the problem of correctly classifying 
aN random input-output associations using a single-layer network, or a tree-like two-layer network, with 
N synapses, up to a value of a close to the theoretical upper bound. For the single-layer case, the 
theoretical bound for binary synapses is Oc — 0.83 [12], while the algorithmic bound as estimated from 
extensive simulations up to = 10® is asp — 0.74. Two more algorithms, obtained as crudely simplified 
versions of the reinforced BP, were later shown dSlIIZI to be able to achieve very similar performances, 
despite being simpler and working in an on-line fashion. The time complexity of all these algorithms was 
measured to be of order N \/log N per pattern; the BP algorithm in particular achieves this performance 
thanks to a Gaussian approximation which is valid at large N. 

When considering multi-layer networks, the original BP approach of [2] can only effectively deal with 
tree-like network structures; fully-connected structures (such as those commonly used in machine learning 
tasks) can not be addressed (at least not straightforwardly) with this approach, due to strong correlations 
arising from a permutation symmetry which emerges in the second layer. 

In this paper, we present a new algorithm for addressing the problem of supervised training of network 
with binary synapses. The algorithm is a variant of the so-called Max-Sum algorithm (MS) [T5| with an 
additional reinforcement term (analogous to the reinforcement term used in |14j). The MS algorithm is 
a particular zero-temperature limit of the BP algorithm; but it should be noted that this limit can be 
taken in different ways. In particular, the BP approach in m was applied directly at zero temperature 
as patterns had to be learned with no errors. In the MS approach we present here, in addition to hard 
constraints imposing that no errors are made on the training set, we add external fields with a temperature 
that goes to zero in a second step. Random small external fields also break the permutation symmetry 
for multi-layer networks. 

In the MS approach, the Gaussian approximation which is used in the BP approach can not be used, 
and a full convolution needs to be computed instead; this in principle would add a factor of N'^ to the 
time complexity, but, as we shall show, the exploitation of the convexity properties of the problem allows 
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to simplify this computation, reducing the additional factor to just log A^. 

This reinforced MS algorithm has very similar performance to the reinforced BP algorithm on single 
layer networks in terms of storage capacity and of time required to perform each node update; however, 
the number of updates required to reach a solution scales polynomially with N, thus degrading the overall 
scaling. On fully-connected multi-layer networks, the MS algorithm performs noticeably better than BP. 

The rest of the paper is organized as follows: in Section |II] we present the network model and the 


mathematical problem of learning. In Section III we present the MS approach for discrete weights. We 
show how the inherent equations can be solved efficiently thanks to properties of the convolution of concave 
piecewise-linear functions, and describe in complete detail the implementation for binary weights. Finally, 


in Section IV we show simulation results for the single and two-layer case. 


II. THE NETWORK MODEL 


We consider networks composed of one or more elementary “building blocks” (units), each one having 
a number of discrete weights and a binary transfer function, which classify binary input vectors. Units 
can be arranged as a composed function (in which the output from some units is considered as the input 
of others) in various ways (also called architectures) that is able to produce a classification output from 
each input vector. 

We denote the input vectors as ^ G {“1; +1}^ (where fi is a pattern index) and the 

weights as ^ (where A: is a unit index). In the following, the are assumed to take q 

evenly spaced values; we will then explicitly consider the cases q = 2 with Wj^ G {—1,1} and q = 3 with 
G {—1) 0,1}. The output of the unit is given by: 


= sign 



( 1 ) 


with the convention that sign (0) = 1. 

We will consider two cases: a single layer network and two-layer comitee machine. In single-layer 
networks, also called perceptrons there is a single unit, and therefore we will omit the index k. Fully 
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connected two-layer networks consist of K units in the second layer, each of which receives the same input 
vector and the output function of the device is 



This kind of architecture is also called a committee or consensus machine |19j . When K = 1, this 
reduces to the perceptron case. In a tree-like committee machine the input vectors would not be shared 
among the units; rather, each unit would only have access to a subset of the input vectors, without 
overlap between the units. For a given N, the tree-like architectures are generally less powerful (in terms 
of computational capabilities or storage capacity) than the fully-connected ones, but are easier to train 
m- Intermediate situations between these two extremes are also possible. In fully-connected committee 
machines there is a permutation symmetry in the indices k, since any two machines which only differ by 
a permutation of the second layer’s indices will produce the same output. 

Throughout this paper we will consider supervised contexts, in which each pattern fj. has an associated 
desired output 

In classification (or storage) problems, M = aNK association pairs of input vectors and corre¬ 
sponding desired outputs are extracted from some probability distribution, and the goal is to find a 
set of weights such that V// G {1,... ,aNK} : 

In random generalization problems, the input patterns are still extracted from some probability 
distribution, but the desired outputs are computed from some rule, usually from a teacher device 
{teacher-student problem). The goal then is to learn the rule itself, i.e. to achieve the lowest possible 
error rate when presented with a pattern which was never seen during the training phase. If the teacher’s 
architecture is identical to that of the student device, this can be achieved when the student’s weights 
match those of the teacher (up to a permutation of the units’ indices in the fully-connected case). 

In the following, we will always address the problem of minimizing the error function on the training 
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patterns: 

aN 

^=1 i 

olN / / K / N \\\ 

= E e E sig“ E - E r? (»'‘) 

/2=1 \ Vfc=l \i=l ) ) ) i 

where 0 (x) is the Heaviside step function 0 (x) = 1 if x > 0 and 0 otherwise. The term T^ {^t) 
the role of an external field, and can be used e.g. to implement a regularization scheme; in the following, 
we will always assume it to be concave. For example we can implement li regularization by setting 
T^(Wj^) = —X\Wj^\ where A > 0 is a parameter. The first term of expression Q therefore counts the 
number of misclassified patterns and the second one favours sparser solutions. 

Throughout the paper, all random binary variables are assumed to be extracted from an unbiased 
i.i.d. distribution. 

Under these conditions, it is known that in the limit of ^ 1 there are phase transitions at particular 
values of a. For single units (perceptrons) with binary ±1 synapses, for the classihcation problem, the 
minimum number of errors is typically 0 up to ac — 0.83. For the generalization problem, the number of 
devices which are compatible with the training set is larger than 1 up to axs — 1.245, after which the 
teacher perceptron becomes the only solution to the problem. 

III. THE MAX-SUM ALGORITHM 

Following |14] . we can represent the optimization problem of finding the zeros of the first term of 
eq. on a complete bipartite factor graph. Starting from the single-layer case, the graph has N vertices 
(variable nodes) representing the Wi values and aN factor nodes representing the error terms (W). 

The standard MS equations for this graph involve two kind of messages associated with each edge of 
the graph; we indicate with the message directed from node /r to variable i at time step t, and 

with (Wj) the message directed in the opposite direction. 

These messages represent a certain zero-temperature limit of BP messages, but have also a direct inter¬ 
pretation in terms of energy shifts of modified systems. Disregarding an insubstantial additive constant. 


message (Wi) represents the negated energy Q restricted to solutions taking a specific value Wi 

for variable i, on a modified system in which the energy depends on Wi only through the factor node 
/i, i.e. in which all terms for all 7 ^ /r are removed from the energy expression (§. Similarly, 

message (kk*) represents an analogous negated energy on a modified system in which the term 

is removed. For factor graphs that are acyclic, the MS equations can be thought of as the dynamic 
programming exact solution. In our case, the factor graph, being complete bipartite, is far from being 
acyclic and the equations are only approximate. For BP, the approximation is equivalent to the one of 
the Thouless-Anderson-Palmers equations | 21 ] and is expected to be exact in the single-layer case below 
the critical capacity m- For a complete description of the MS equations and algorithm, see m- 
The MS equations m for energy (§ are: 




max 

{Wj}.^.E^iW)=0 





'7t + l 


( 3 ) 

( 4 ) 


where and Zl^^ are normalization scalars that ensure Yhwi (kki) ~ Ewi i^i) — 0 ^''^4 

can be computed after the rest of the RHS. At any given time t, we can compute the single-site quantities 


(Wi) = Ti (Wi) + (W,) - Z( (5) 

and use them to produce an assignment of the W’s: 


IF/ = argmax^y^T* (IF*) ( 6 ) 

The standard MS procedure thus consists in initializing the messages (kF*), iterating eqs. Q and 
@ and, at each time step t, computing a vector IF* according to eqs. and Q until either E (IF*) = 0 
(in the absence of prior terms, i.e. when A = 0), or the messages converge to a fixed point, or some 
maximum iteration limit is reached. 

Strictly speaking, standard MS is only guaranteed to reach a fixed point if the factor graph is acyclic, 
which is clearly not the case here. Furthermore, if the problem has more than one solution (ground state). 
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the assignment in eq. Q would not yield a solution even in the acyclic case. In order to (heuristically) 
overcome these problems, we add a time-dependent reinforcement term to eqs. Q and ([^, analogously 
to what is done for BP |14j : 






{W,)=rt m) + r, (IP,) + E 


t 




4/* (IP,) = (IP,) + r, (IP,) + (IPO - z* 


( 7 ) 

( 8 ) 


where r > 0 controls the reinforcement speed. This reinforcement term in the case of standard BP 
implements a sort of “soft decimation” process, in which single variable marginals are iteratively pushed 
to concentrate on a single value. For the case of MS, this process is useful to aid convergence: on a system 
in which the MS equations do not converge, the computed MS local fields still give some information about 
the ground states and can be used to iteratively “bootstrap” the system into one with very large external 
helds, i.e. fully polarized on a single configuration |22) . The addition of this term introduces a dependence 
on the initial condition. Experimentally, by reducing r this dependence can be made arbitrarily small. 


leading to more accurate results (see Sec. IV), at the cost of increasing the number of steps required for 
convergence; our tests show that the convergence time scales as r~^. 

Furthermore, in order to break symmetries between competing configurations, we add a small 
symmetry-breaking concave noise P) (IP,) <C 1 to the external fields Tj (IP,); this, together with the 
addition of the reinforcement term, is sufficient to ensure — for all practical purposes — that the argmax 
in ([^ is unique at every step of the iteration. 


A. Max Convolution 

While Eq. 0 can be efficiently computed in a straightforward way, the first term of Eq. 0 involves 
a maximum over an exponentially large set. The computation of Eq. Q can be rendered tractable by 
adding 0 = maxA L ^A, where L (x, y) is equal to 0 if x = y and —oo otherwise, which leads 

to the following transformations; 



10 


I't"-!. (Wi) + z‘+‘, = 




max 




E '^5^,. ("'i 


max < maxL A, 


max < max > '^',-s./, (Wj 

A:a^(A+CWi)>0 ] i^Wj=A ^ 


max j (A), 

A:a^(A+^fVKi)>0 


(9) 


where in the last step above (A) is dehned as: 


(A) = “ax I . 




( 10 ) 


The right-hand side of (10) is usually called a “Max-Convolution” of the functions fj{Sj) = 


foi" j / h is analogous to the standard convolution but with operations (max,-)-) 
substituting the usual (-I-, x). As standard convolution, the operation is associative, which allows to 
compute the convolution of the N — 1 functions in a recursive way. As the convolution of two functions 
with discrete domains {0,..., qi} and {0,..., ( 72 } respectively can be computed in qiq 2 operations and 
has domain in {0,..., gi -|- q 2 }, it follows that (10) can be computed in O (A^) operations. In principle, 
this computation must be performed N times for each pattern /r to compute all messages, in a total 
of time O (A^). 

A technique like the one described in [23H25] can be employed to reduce this by a factor A, coming 
back again to O (A^) operations per pattern update, as follows. Precomputing the partial convolutions 
Ln and Rn of /i,..., /„ and fn, ■■■, In (respectively) for every n in 1,..., A can be done using A^ 
operations in total; then the convolution of can be computed as the convolution of Lj-i and Ri+i. 

Computing this convolution would require 0{N‘^) operations but fortunately this will not be needed; the 
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computation of Q can proceed as: 

(Wi) + = max I max L^.i (z) + Ri+i (A - z) 

A:a%{A+^I^Wi)>0 2 

= max < Li-i (z) + max Ri+i (A — z) 

^ [ A-.a'^{A+i'^Wi)>0 

= max |Lj_i (z) + {z + ( 11 ) 

where we defined R^ (x) = max^.o-(^+a,)>o-R* (A). As the vectors R^ can be pre-computed recursively 
in a total time of O {N‘^) and © requires time O (N), we obtain a grand total O (-/V^) operations per 
pattern update, or 0{MN‘^) per iteration. Unfortunately, this scaling is normally still too slow to be of 
practical use for moderately large values of N and we will thus proceed differently by exploiting convexity 
properties. However, note that the above scaling is still the best we can achieve for the general case in 
which regularization terms are not concave. 

At variance with standard discrete convolution, in general Max-Convolution does not have an analogous 
to the Fast Fourier Transform, that would allow a reduction of the computation time of a convolution of 
functions with 0{N) values from to A log A. Nevertheless, for concave functions the convolution can 
be computed efficiently, as we will show below. Note that for this class of functions, an operation that is 
analogous to the Fast Fourier Transform is the Legendre-Fenchel transform [26j , though it will be simpler 
to work with the convolution directly in the original function space. 

First, let us recall well-known results about max-convolution of concave piecewise-linear functions in 
the family C = {/ : M>o {—oo}} |26j. First, the max-convolution /12 '^= fi & f 2 of /i ,/2 £ C 

belongs to C. Moreover, /i ® /2 can be built in an efficient way from /i and f 2 - Start with '^= 
inf {x : (/i ® 72 ) (x) > — 00 }, which is easily computed as = x\ + x\ with x\ = inf {x : fi (x) > — 00 }. 

Moreover, (/i ® 72 ) (x^) = 7i (®i) + 72 (^i)- Then, order the set of linear pieces from 7i and 72 in 
decreasing order of slope and place them in order, starting from (x^, f (x^)) to form a piecewise-linear 
continuous function. The method is sketched in Fig. In symbols, let us write each concave piecewise- 
linear function fi (x) when x > x^ as: 

hi 

fi (x) = ^(x- + fi (^l) * = 1, 2 
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with a®- G [0, oo] for j = 2,..., fc* and x\ < X 2 < ■ ■ ■ < x\,. Here we used the notation[27j ^ (|?/| + y). 

This function is concave, as for x € Xj,Xj_^_i the slope is 6* = Ylk=i^k clearly decreasing with 

j. To compute the convolution of /i and / 2 , just order the slopes by, i.e. compute a one to one map 
TT : {i,j) I—>■ c from couples i G {1, 2},1 < j < ki to integers 1 < c < ki + k 2 such that tt {i,j) < vr {i',j') 
implies 6®- > 6*/. The max convolution for x > x^ is still concave and piecewise-linear, and thus it can 
be written as: 


ki2 

(/i ® /2) (x) = X] + (/i ® /2) (xP) 

C=1 

where ki 2 = ki + k 2 — 1- For each c we can retrieve (z (c) ,j (c)) = 7r~^ (c); with this, the parameters of 
the convolution are = 6®.W, and x^l^ = x^ + 

For more details about the max-convolution of piecewise-linear concave functions, see e.g. |26l Part 
II]. 

We now consider the case of functions defined on a discrete domain. Let /, g be concave discrete 
functions in 


F* = {/ : {0,..., g - 1} M U {- 00 }} 

We will define the continuous extension / G C as the piecewise-linear interpolation of /, with value —00 
for arguments in (—oo,0) L) {q — 1 , 00 ). This can be written as: 

<? 

/ (a:) = ^{x-j + 1 )+ aj + f ( 0 ) 
i=i 

with ai = / (1) - / (0), aj+i = f {j + 1) - 2/ (j) -h f{j - 1) (implying Ug = - 00 ). It is easy to see that 
h = f & g coincides with the discrete convolution of / and g in its (discrete) domain; the reason is simply 
that h is also piecewise-linear, with kinks only in discrete values ,{0,... 2 (g — I)}. 

When computing the convolution of N functions /i,..., /^r with domain in {0,..., g — I}, one can 
clearly order all slopes together in one initial pass, using Nq log {Nq) comparisons. It is also easy to see 
that, if one has the full convolution, it is simple to compute a “cavity” convolution in which one function /j 
is omitted, in O {q) time: this is achieved by simply removing the q slopes of /* from the full convolution. 
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FIG. 1. Sketch of the discrete max-convolution f&goi piecewise-linear concave functions /, g. The result is simply 
obtained by sorting the pieces of the two functions in descending order of slope. 


In order to apply this to eq. (10) the only remaining step is a trivial afhne mapping of the functions 
arguments on a domain {A + tB : t £ {0,..., q — 1}} to the domain {0,..., g — 1}. In the following, 
we will show explicitly how to do this for the binary case q = 2, but the argument can be easily generalized 
to arbitrary q. Note that, while in the binary case the 'h functions are linear and thus trivially concave, in 
the general case we need to ensure that both the initial values (Wj) and the external fields Fj (Wj) 
are concave; in such case, the iteration equations Q, Q and Q ensure that the concavity property holds 
for all time steps t > 0. 


B. The binary case 


We will show explicitly how to perform efficiently the computations for the binary case. In this 
case we can simplify the notation by parametrizing the message with a single scalar, i.e. we can write 




and (Wi) = (^*) = Eqs. @ and Q then become: 




- ( max Wjip 

2\w-.Wi=lAEJW)=0^ 


t 


max > 

W:Wi=-lAEJW)=0^ 




( 12 ) 




t-i 


+ 7 .+ E"* 






(13) 
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Correspondingly, eqs. 


and Q simplify to: 

i’i=rt + li + '^4>\ 
wl = signV’- 




(14) 

(15) 


In order to apply the results of the previous sections, and perform efficiently the trace over all possible 


assignments of W of eq. (12), we first introduce the auxiliary quantities 




Ti (A) = max F* (S) 
{S:E,5,=a} 


(16) 

(17) 


For simplicity of notation, we will temporarily drop the indices \x and t. We will also assume that 
all values are different: as remarked above, the presence of term Fj is sufficient to ensure that 

this is the case, and otherwise we can impose an arbitrary order without loss of generality. With this 
assumption, the function F, which is defined over A = {—N, —N + 2,..., N — 2, N}, has a single absolute 
maximum, and is indeed concave. The absolute maximum is obtained with the special configuration 
S = argmax| 5 j F (S), which is trivially obtained by setting Si = for all i. This configuration 

corresponds to a value A = variable flip with respect to this conhguration, i.e. any i for 


which Si = —Si, adds a “cost” AF,- = 2 






i in two groups 5+ and 5_ defined by = ■ 


in terms of F (S). Therefore, if we partition the indices 
i : .Sj = ±l|, and we sort the indices within each group 


in ascending order according to AF,, we can compute the function F (A) for each A by constructively 
computing the corresponding optimal configuration S, in the following way: we start from F then 
proceed in steps of 2 in both directions subtracting the values AFj in ascending order, using the variable 
indices in S^ for A < A and those in S- for A > A. 

This procedure also associates a “turning point” Ti to each index i, defined as the value of A for which 
the optimal value of Wi changes sign, or equivalently such that F ^Tj + Si^ — F — Si^ = AFj. This 
also implies that: 


F (A) = F (a) - ^ 0 (F - A)) AF 


(18) 
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We can also bijectively associate an index to each value of Tj, by defining jk such that Tj^ = k. 
Next, consider the same quantity where a variable i is left out of the sum (see eq. ®) 

(S«) = 

j¥=i 

J-® (A) = 




max 


(19) 

( 20 ) 


Clearly, one gets the same overall picture as before, except with a shifted argmax, and shifted turning 
points. The shifts can be easily expressed in terms of the previous quantities, and the expressions used 


for computing eq. (12) as: 




2 


max 

A;(T^A>0 


jr(0 (A - 1) - 


max 
A: cr^ A>0 


(A + 1) 


( 21 ) 


The full details of the computation are provided in the Appendix, Sec. |Vj Here, we report the end 
result; 


= ii (© {^d) 0 (-A + s, + l) (0 m - 1) h,, + 0 i-Ti + 1) h,,) + (22) 

+ 0 (-CT^) 0 (a - + l) (0 (Ti + 1) hj_, + 0 (-T, - 1) 

where hj = From this expression, we see that we can update the cavity fields very 

efficiently for all i, using the following procedure: 


• We do one pass of the whole array of hi by which we determine the Si values, we split the indices 
j into 5"+ and S- and we compute A. This requires O (N) operations (all of which are trivial). 


• We separately partially sort the indices in 5*+ and S- and get j_ 2 , jo and j 2 and the turning points 
Ti. This requires at most O (NlogN) operations. Note that we can use a partial sort because we 
computed A, and so we know how many indices we need to sort, and from which set S±, until we 
get to the ones with turning points around 0; also, we are only interested in computing 0 (Tj — 1) 
and 0 {Ti + 1) instead of all values of Ti. This makes it likely for the procedure to be significantly 
less computationally expensive than the worst case scenario. 

• For each i we compute (j)^^i from the equation above. This requires 0(1) operations (implemented 
in practice with three conditionals and a lookup). 
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IV. NUMERICAL RESULTS 

We tested extensively the binary case q = 2 with Wi G {—1,+1} and the ternary case q = 3 with 
Wi G { — 1,0,1}, for single layer networks. 

We start from the binary case. Fig. shows the probability of finding a solution when fixing the 
reinforcement rate r, for different values of r and a. Reducing r allows to reach higher values of a; the 
shape of the curves suggest that in the limit r —)• 0 there would be sharp transitions at critical values of 
a’s. In the classification case. Fig. [^, the transition is around a ~ 0.75, while the theoretical critical 
capacity is ac = 0.83. This value is comparable to the one obtained with the reinforced BP algorithm 
of [TTj. In the generalization case, there are two transitions: the first one occurs around a ~ 1.1, before 
the first-order transition at aTS = 1.245 where, according to the theory, the only solution is the teacher; 
the second transition occurs around a ~ 1.5. This second transition is compatible with the end of the 
meta-stable regime (see e.g. my, indeed, after this point the algorithm is able to correctly infer the teacher 
perceptron. 

A second batch of experiments on the same architecture, in the classification case, is shown in Fig. 

In this case, we estimated the maximum value of r which allows to find a solution, at different values of 
N and a; i.e. for each test sample we started from a high value of r (e.g. r = 10“^) and checked if the 
algorithm was able to find a solution; if the algorithm failed, we reduced r and tried the same sample 
again. In the cases shown, the solution was always found eventually. The results indicate that the value 
of r required decreases with N, and the behaviour is well described by a power low, i.e. r = aN^ with 
a < 0 and 6 < 0, where the values of a and b depend on a. Since the number of iterations required is 
inversely proportional to r (not shown), this implies that the overall solving time of the MS algorithm is 
of O log (N)^, i.e. it is worse than the reinforced BP in this respect. The value of b is between 0 

and —0.5 up to a = 0.6, after which its magnitude decreases abruptly (see Fig. [^). The behaviour for 
large a seems to be reasonably well fit by a curve b(a) = , suggesting the presence of a vertical 

asymptote at au = 0.755 ± 0.004, which is an estimate of the critical capacity of the algorithm in the 
limit of large N. 

In the ternary single layer case, we tested learning of {—1,1} random patterns with ternary {—1, 0,1} 
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FIG. 2. Solving probability. Probability of finding a solution for different values of a, in the binary perceptron 
case with N = 1001, with different values of the reinforcement rate parameter r. Performance improves with lower 
values of r. A. Classification case, 100 samples per point. The theoretical capacity is Oc = 0.83 in this case. B. 
Generalization case, 20 samples per point. In this case, the problem has multiple solutions up to ars = 1.245, after 
which the only solution is the teacher. 


weights and concave bias (i.e. prior). In practice, we use the function Fj (Wj) = X5 {Wi) + r' (Wj) (where 
r' is the symmetry-breaking noise term and A is sufficiently large) to favour zero weights, so solutions 
with a minimimal number of zeros are searched, i.e. we add an Iq regularization term. Results (See 
Fig.U) are qualitatively similar to the {—1,1} case with a larger capacity (around a = 1; the critical 
capacity is Oc = 1-17 in this case). The average non-zero weights in a solution grows when getting closer 
to the critical a up to a value that is smaller than 2/3 (the value that makes the entropy of unconstrained 
{—1,0,1} perceptrons largest). 

In the fully-connected multi-layer case, the algorithm does not get as close to the critical capacity 
as for the single-layer case, but it is still able to achieve non-zero capacity in rather large instances. 
For example, in the classification case with binary synapses, N = 1001 inputs, K = 3 hidden units, 
the algorithmic critical capacity is a ~ 0.33 when r = 10“® (tested on 20 samples), corresponding to 
storing M = 1001 patterns with 3003 weights (thus demonstrating a greater discriminatory power than 
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FIG. 3. Maximum reinforcement rate. A. Average value of the maximum reinforcement rate r which allows 
to find a solution, in the binary perceptron classification case, at various values of N and a, in log-log scale. The 
reinforcement rate decreases with N and a. Error bars show the standard deviation of the distributions. Black lines 
show the result of fits of the form r {a, N) = a (a) one for each value of a. The number of samples varied 

between 100 for N = 1001 and 10 for N = 32001. B. The fitted values of the exponents b{a) in panel A. The 
continuous curve shows a fit of the data for a > 0.6 by the function b{a) = ■ The fit yields c = —0.063 ±0.002 

and au = 0.755 ± 0.004. The value of ajj is an estimate of the critical capacity of the algorithm. 


the single-layer case with the same input size). The reason for the increased difficulty in this case is not 
completely clear: we speculate that it is both due to the permutation symmetry between the hidden units 
and to replica-symmetry-breaking effects: these effects tend to trap the algorithm — in its intermediate 
steps — in states which mix different clusters of solutions, making convergence difficult. Still, the use 
of symmetry-breaking noise helps achieving non trivial results even in this case, which constitutes an 
improvement with respect to the standard BP algorithm. 


V. CONCLUSIONS 


Up to now, the large N limit could be exploited on BP equations for the learning problem with dis¬ 
crete synapses to obtain an extremely simple set of approximated equations that made the computation 
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FIG. 4. Learning of random {—1,1} patterns with a ternary w = {—1,0,1} perceptron, with dilution (regularization) 
prior term; N = 1001, 20 samples per point. For solved instances with r = 10“^^, the average fraction of non-zero 
weights is also shown (standard deviations smaller than point size). 


of an iteration to scale linearly with the problem size MN. For the MS equations however, those ap¬ 
proximations cannot be made and a naive approximation scales as MN^ which is normally too slow for 
most practical purposes. In this work, we showed that MS equations can be computed exactly (with no 
approximations) in time MNlogN, rendering the approach extremely interesting in practice. A word 
is in order about the MS equations with reinforcement term, which we propose as a valid alternative to 
Belief Propagations-based methods. Although we cannot claim any universal property of the reinforced 
equations from theoretical arguments and we only tested a limited number of cases, extensive simula¬ 
tions for these cases and previous results obtained by applying the same technique to other optimization 
problems of very different nature [2211^5113Uj have confirmed the same qualitative behaviour; that is, that 
the number of iterations until convergence scales as r~^ and that results monotonically improve as r 
decreases. As an additional advantage of MS, inherent symmetries present in the original system are nat¬ 
urally broken thanks to ad-hoc noise terms that are standard in MS. The MS equations are additionally 
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computationally simpler because they normally require only sum and max operations, in contrast with 
hyperbolic trigonometric functions required by BP equations. Extensive simulations for discrete {—1,1} 
and {—1,0,1} weights show that the performance is indeed very good, and the algorithm achieves a 
capacity close to the theoretical one (very similar to the one of Belief Propagation). 
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APPENDIX 


Details of the computation of the cavity fields. 


In this section we provide the full details of the computation leading to eq. (22). 

As noted in the main text, the expression of the cavity quantities (A) (see eq. (20)) is analogous 
to that of the non-cavity counterpart T (A) (eq. where the argmax has changed to A^*^ = A — Si, 

and the turning points have changed: 


J-W (A) = .fW (^aW) - ^ 0 (-Sj - a)) AFj (23) 

We need to express the relationship between the old turning points and the new ones: having omitted 
the variable i, it means that there is a global shift of —Si, and that the turning points to the left (right) 
of Ti have shifted to the right (left) if 5* = 1 (5* = — 1, respectively): 

T® = Tj - Si + 2SiQ (^Si {Ti - Tj)^ 

= Tj+sign (Ti-Tj-S^'j (24) 

(note that we chose to use the convention that 0 (0) = 0). 

Therefore we obtain: 


^(^) = ^(i) ('aW) - ^ 0 (4- (r,- + sign (Ti - Tj - Si') - a)) AFj 


(25) 
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Next, we consider the cavity quantity: 


Ci (A, Si) = max V 

= E(') (A - S^) 


(26) 


which allows us to write eq. (12) as 




^ L - . max Ci (A, -4^) 


, A;(T^A>0 


A:o-^A>0 


S2 

2 


max 

A;(T^A>0 


jr(d (A - 1) - 


max 
A: cr^ A>0 


(A + 1) 


(27) 


(this is eq. (21) in the main text). 

Note that (A) is concave and has a maximum at A^*) = A — Sj. Using this fact, and eq. (25), we 
can derive explicit formulas for the expressions which appear in the cavity field, by considering the two 
cases for separately, and simplifying the result with simple algebraic manipulations afterwards: 


^ m^ax^^E^*) (A - 1 ) = 0 (a^) (o (a - + l) (a - 5 *) + 0 (-A + 5 * - l) ( 0 )) + 

+0 (-a^) (0 (-A + Si- 1 ^ (a - + 0 (a - + 1) (-2)) 

= 0 (a - + 1)) e(') (a - + 

+0 (cT^) 0 (-A + 5,-1) ( 0 ) + 

+0 (-(T^) 0 (a - 5 , + 1) (-2) 


max (A + 1) 
A:a^A>0 


0 (a^) (0 (a - 5 , - 1) (a - 5 *) + 0 (-A + + l) (2)) + 

+0 (-u^) (0 (-A + + 1) e(') (a - 5,) + 0 (a - - 1) e(') ( 0 )) 

0 (a - - 1)) (a - + 

+0 (u^) 0 (-A + + 1) ( 2 ) + 

+0 (-u^) 0 (a - - 1) E(') (0) 
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Plugging these back in the expression for the cavity held, we can reach — again by simple algebraic 
manipulations — an expression which only uses (~2), (0) and (2): 

Ci = I ((® (<"& (a - s* + l)) - ® (^ - -S. - 0)) - •*■) + 

+0 (0 (-A + Si- 1 ^ ( 0 ) - 0 (-A + 5i + 1 ) ( 2 )) + 

+ 0 (-a^) (0 (a - + 1 ) J-(*) (- 2 ) - 0 (^A - - 1 ) ( 0 )) ) 

= ^ (5(A,5,)sign(a^)^«(0) + 

+0 (a^) (-,5 (a, (2 ) + 0 (-A + 5^ - l) (0) - (2))) + 

+ 0 {-P^) (<5 (a, 5*) P^^ (- 2 ) + 0 (a - 5i - 1 ) (- 2 ) - P^^ (0)) ) ) 

= ^ [5 (A, (sign (a^) (^(*) (0) - (2 ct^))) + 

+0 {P^) 0 (-A + Wip - 1 ) (0) - J-W (2)) + 

+ 0 i-P^) 0 (A - mp - 1 ) (- 2 ) - ( 0 ))) 

= y (0 (^d) 0 (-^ + Si + 1^ (^P^'> (0) - P^'> (2)) + 

+ 0 {-P^) 0 (A - + 1 ) (- 2 ) - ( 0 ))) 

These expressions can be further simplihed, since the differences between the values of P^'^ at neigh¬ 
boring values only depends on the “steps” induced by the spins which are associated with turning points 
in that region: 


(0) - pi) (2) = (0 ( 5 ,- (y + sign (Ti - Tj - 5,) - 2 )) - 0 (y (y + sign (t, - T, - 5,)))) Ay 

= - E 1 - ^*)) 

= -4 (1 - <5 (i, jo)} 0 (y - y) Ay„ - 4 (1 - <5 (i, J 2 )) 0 (-F + 2 + y ) AFj, 

= -40 m - 1) Ay-o - 40 (-y +1) Ay, 
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where in the last step we used the Kronecker deltas to get rid for the differences between the two cases 
for Si- The other case is very similar: 


jr(i) (_2) _ ( 0 ) = (0 ( 5 ,- (r,- + sign (t, - T, - 5,))) - 0 (^Sj (Tj + sign (t, - T,- - 5,) + 2 ))) AFj 

j¥=i 

= {Tj, -1 - sign (t, - T, - 5.)) AF, 

= -Sj_, (1 - 5 e(Ti + 2- Si) AF,_, - 4 (1 - 5 (i, jo)) 0 [Si - T,) AFj, 

= -4,0 m + 1) - 40 (-1 - T,) AF ,, 


Going back to the cavity fields, and defining hj = —kSjAFj = —^44 , we finally get eq. (22): 


44* - 


= ef (0 (4) 0 (-A + s, + i) (0 m - 1) 4 + 0 i-Ti +1) 4) + 

+ 0 (-4) e(^A-Si + l) (0 {Ti + 1) hj_, + 0 i-Ti - 1) 4)) 
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